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We investigate a class of stochastic growth models involving com- 
petition between two phases in which one of the phases has a com- 
petitive advantage. The equilibrium populations of the competing 
phases are calculated using a mean field analysis. Regression prob- 
abilities for the extinction of the advantaged phase are calculated 
in a leading order approximation. The results of the calculations 
are in good agreement with simulations carried out on a square 
lattice with periodic boundaries. 

The class of models are variants of the Williams- Bjerknes model 
for the growth of tumours in the basal layer of an epithelium. In 
the limit in which only one of the phases is unstable the class of 
models reduces to the well known variants of the Eden model. 



1 Introduction 



Stochastic pattern formation is ubiquitous in science and technology (see, e.g., 
[1-4] and the many references therein) . Some examples include viscous finger- 
ing in fluids, dendritic crystallization, electro-deposition, dielectric breakdown 
and the growth of bacterial colonies. It is remarkable that, despite their dis- 
parate physical origins, many such patterns can essentially be described by one 
of two generic models. These models are distinguished by the nature of the 
stochastic growth at the interface, which is usually taken to be either: (i) uni- 
form - as in the Eden model [5,6], or (ii) Laplacian - as in the diffusion-limited 
aggregation model [7]. 

These generic models and their numerous refinements have been intensively 
studied and have proved highly successful in describing the growth of patterns 
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in which one stable phase propagates into a second unstable phase. For exam- 
ple, in the formation of snowflakes the snow crystal (stable phase) propagates 
into water vapour (unstable phase). In viscous fingering a low viscosity fluid 
(stable phase) is pumped under pressure into a background high viscosity fluid 
(unstable phase). These are examples of Laplacian growth processes. An ex- 
ample of a uniform growth process is the spread of a bacterial colony (stable 
phase) via cell-division into a nutrient-rich environment (unstable phase) [8]. 

However, many important growth processes have two or more competing un- 
stable phases. An example is the Williams-Bjerknes (WB) model [9] for the 
growth of tumours, in which both the cancerous cells and the normal cells 
are unstable (both may be displaced by cell division). The limiting growth is 
governed by this competition. The cancerous cells have the competitive advan- 
tage of dividing faster than normal cells. Importantly, early simulations of the 
model [9] revealed that the invasive nature of the tumours can be accounted 
for solely by an almost evenly balanced contest between the cancerous cells 
and the normal cells. These early simulations also investigated the regression 
probability, defined as the probability that the advantaged phase, starting out 
as a single seed, will become extinct under the growth process. Subsequent 
simulations [10] and rigorous mathematical studies [11,12] revealed that the 
abnormal region, whenever it survives, can be described as an asymptotic 
shape (with a one-dimensional boundary), the average radius of which grows 
linearly in time. 

In this paper we introduce a class of uniform stochastic growth models (which 
includes the WB model) involving competition between two phases. The equi- 
librium populations of the competing phases are calculated using a mean field 
analysis similar to that used recently for a dynamical model of virus spread 
[13]. Leading order approximations for regression probabilities are also calcu- 
lated. 

The mean field results and regression results are in good agreement with sim- 
ulations carried out on a square lattice with periodic boundaries. 



2 Competing Growth Models - Mean-Field Equations 

Consider a regular lattice with N sites where N A are of type A and N B are type 
B. Define an interface site as one which has one or more nearest neighbouring 
sites of the opposite type. Let Ma an d Mb denote the number of type A 
interface sites and type B interface sites respectively. The total number of 
nearest neighbour sites surrounding a lattice site is the lattice co-ordination 
number denoted here by n c (for example, n c = 4 on the square lattice). 
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In the models below the numbers of type A sites and type B sites change in 
time as a consequence of interface interactions which occur when an interface 
site is changed into the opposite type. Let Pi^ B denote the probability for site 
% to be changed to a type B site and similarly for p^ A . The complete set of 
transition probabilities for the site % are then; 



Pi(B)->A=Pi^ASi,B, (1) 

Pi(A)^B=Pi^B$i,A> ( 2 ) 

Pi(A)^A=Pi->A&i,A, (3) 

Pi(B)->B=Pi^B$i,B- (4) 



where i(A) and i(B) is used to represent the situation in which site % is a type 
A site and a type B site respectively. The mean field probability that a B site 
will be changed into an A site is thus 

N 

PB^A = ^2pi^A5i,B- (5) 



Similarly 

N 

PA->B = J2Pi^B<)i,A (6) 

i 



is the mean field probability that an A site will be changed into a B site. Using 
these probabilities we can write down general mean field population equations 
for this class of models which relates the numbers of species before (time t) 
and after (time t + 1) an interface interaction: 



N A {t + 1) =p B ^ A {t) - p A ^ B (t) + N A (t), (7) 
N B (t + 1) =p A ^ B {t) - p B ^ A {t) + N B (t). (8) 

The mean field steady state defined by N A (t + 1) = N A (t) = N A and N B (t + 
1) = N B (t) = N B thus yields the (equivalent) equilibrium conditions 

P* B -> A = P* A ^ B , (9) 

N N 

EpU^ = Ep^- (io) 

i i 

We shall make use of both forms of the mean-field equilibrium condition in 
the following. Note that the steady state mean field probabilities p* B ^ A and 
p* A ^ B are in general functions of N A ,Af A ,N B ,Af B . 
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3 Williams-Bjerknes Model 



In the WB model, the spread of cancers is modelled by the growth of two 
types of cells, cancerous cells and normal cells, which compete via cell division 
in a basal layer. The cancerous cells are taken to divide k times faster than 
normal cells. This factor k is called the carcinogenic advantage. It is supposed 
that when a cell in the basal layer divides the daughter cell displaces a neigh- 
bouring cell up out of the basal layer [14]. The implementation of the model 
is straightforward. Start with a regular lattice of sites (basal layer) labelled 
type A (cancerous) or type B (normal). Select a site at random with a bias 
(due to the carcinogenic advantage) for selecting type A. Choose one of the 
neighbouring sites of the selected site at random and convert it to the same 
type as the selected site. Typically there are many growth steps that do not 
change the configuration of the basal layer. 

We have simulated this model within a rectangular region covering [0, 100] x 
[—100,100] square lattice sites using periodic boundaries in both directions. 
Initially the top half of the rectangular region is comprised of abnormal cells 
and the bottom half is comprised of the normal cells. The configurations of 
normal and abnormal cells for the case k — 2 are shown at four different time 
snapshots in Fig. 1. Two points of note are: i) the interface between normal 
and abnormal cells is highly irregular; ii) the ratio abnormal cells to normal 
cells increases in time until (on average) the advantaged species completely 
dominates. 



3. 1 Mean Field Equilibrium 



To determine the mean field equilibrium state for this model we first define the 
model transition probabilities. The probability that a site % will be changed 
(B — > A) or reconfigured (A — > A) into type A is 




where rii(A) is the number of neighbouring A sites adjoining site i. The first 
factor is the carcinogenic advantage for a division of an A cell. The second 
factor is the probability that the A cell chosen for division is one of the neigh- 
bours of site i. The final factor is the probability that if a neighbouring A cell 
for site i is chosen for division it will displace the cell at site i out of the basal 
layer leaving an A cell at that location. Similarly the probability that a site i 
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will be changed or reconfigured into type B is 



The normalization of the probabilities, 



N 



P = '£(p i ^ A + p i ^ B ) = l, (13) 



follows immediately from the identities 



Y,^ = N A , (14) 

i n c 
I It. 



Substituting the expressions for the transition probabilities (11), (12) into the 
mean field equilibrium condition (10) yields 

N N 

KY j n l {A)5 hB = Y j n l {B)5 lA . (16) 



This result can be simplified by the identity 

N N 

Y / n t (A)5 i , B = Y / n i (B)5 hA ; (17) 



which follows from the Eqs (14), (15), together with the neighbour summation 
rule rii(A) + rii(B) = n c . Using the identity (17) the mean field equilibrium 
result (16) simplifies to 

k=1. (18) 



From our stochastic simulations described above we have have obtained a plot 
in Fig. 2 of the time taken for all B cells to die out (k > 1) and the time for all 
A cells to die out (k < 1) as a function of log/t. The extinction time rapidly 
increases with decreasing \k — 1| and there is a clear possibility for a steady 
state with both species at k — 1. 
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A different form of the equilibrium condition can be derived by writing the 
mean field probability that any B type cell is changed into an A type cell as 



= (s^nnd (]£) < 19 > 



The first factor on the right is the probability that an A cell divides, the second 
factor is the probability that the dividing A cell is an interface cell and the 
third factor is the probability that A interface cells displace B interface cells 
(there is a finite probability that A cells might divide without any resultant 
change in the configuration). Similarly the mean field probability that any A 
type cell is changed into a B type cell can be written as 



kN a + N b J \N b 



Substituting (19) and (20) into the mean field equilibrium condition (9) now 
yields 

kNaocb^a = Nbola->b- (21) 



Comparing (21) with the previous equilibrium condition (18) we obtain the 
additional equilibrium result 

— = ¥r- (22) 



Note that in general we can write 

N 

Y j n l (B)5 lA = a B M B , (23) 



where the factor as is included because some of the B interface cells are 
counted more than once in the sum, and similarly 

N 

J2ni(A)5 i:B = a A M A - (24) 



Combining (23) and (24) with the identity (17) thus yields 

a A Af A = a B Af B . (25) 
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Comparing the general result (25) with the equilibrium condition (22) we now 
have 

OA = OB^A (26) 
Ct B CtA^B 



3.2 Regression 



Here we consider the case in which a single cancerous site (Na = 1) is sur- 
rounded by normal sites (N B = N — 1). The transition probabilities are thus 



Pi ^ A = L + N-l ) nM ^ B > (27) 

/ N — 1 \ fni(B)\ 1 r 

= (29) 

Summing over the lattice sites we have 



p — Gdbr)' (31) 

'^fa^l)' (32) 

Pa^a = 0, (33) 
' N-2 ' 

.K + N-l, 

The probability Pa^b is the probability for regression in one step. The prob- 
ability for regression without any growth within m steps is given by 



PA^Es(rn) = ^2(p B ^ B ) k 1 p A ^B, (35) 
k=l 

1 _ ( jV=g \ m 



/</; /; ( \. 2 , ) • (31) 



K+ 1 



(36) 



and the probability for regression without growth after infinitely many steps 
is 
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cc 



PA^B(oo) = J2(PB^B) k 1 PA^B, (37) 
k=l 
1 



K+ 1 



(38) 



We anticipate that the probability for regression without growth will be the 
dominant term in the overall probability for regression, which includes the 
probability for regression after growth, so that (38) may be taken as a leading 
order approximation. To investigate this further consider regression within 
three steps. The possible scenarios and associated probabilities are: 



A—>B;B 
B -> B;A 
B -> B;B 
B -> A: A 



B ; B 
B ; B 
B-A 
B ; A 



B 
B 
B 
B 



Pa^b] 
Pb->bPa^b] 
(Pb^b) 2 Pa^b; 
Pb^aPa^bPa^b; 



(39) 
(40) 
(41) 
(42) 



The sum of the probabilities in the first three scenarios above is the probability 
for regression without growth within three steps (Eq (35) with m — 3) and 



Pa^b = I - 



iV-2 



2k + N -2/ \N -2 



(43) 



is the probability A — > B in a cluster with two adjacent type A cells. It 
immediately follows that the ratio the probability for regression after growth 
within three steps to the probability for regression without growth within 
three steps scales as 1/iV 3 . 

We have studied regression in simulations with k — 2 over [—10, 10] x [—10, 10] 
sites on a square lattice. In Fig. 3 we have plotted the frequency of regres- 
sion within three, five and ten steps versus the total number of runs used in 
the frequency estimate, over a range of the total number of runs up to 10 6 . 
The upper two horizontal lines on the plot are the probabilites for regression 
without growth after five and ten steps and the lower horizontal line is the 
full probability for regression after three steps. The comparison between the 
Monte Carlo results and the algebraic regression results provides evidence that 
regression without growth dominates the regression process. 



3.3 Limit k — > oo 



In the limit k — > oo the B cells never divide and the transition probabilities 
become: 
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Pi(A)-+B = 0, 
Pi(B)~*B = 0, 

Pi(A)-^A 



'n t (A)\ 1 
N A n c 



Si,A, 



(m(A)\ l 

Pi(B)^A = | — — ] — <Ji,B- 



n. 



(44) 
(45) 

(46) 
(47) 



In this case there are only two possible events A — > A and B —> A. The event 
A — > A affects the time scale for configurational changes but does not affect 
the configurations themselves. Hence if we do not concern ourselves with the 
time scales for change we can set Pi(A)->A — and renormalize 

ni(A)5 iyA 

Pi(B)^A - (48) 



where the sum is over all type B interface sites. In this limit the model is 
equivalent to the original model introduced by Eden [5,6], which is also referred 
to as the Eden B model [15]. 



4 Interface Model 



In this model only interface sites are considered for division. Select an inter- 
face site at random with a bias for selecting type A and convert one of the 
neighbouring sites of the selected site to a site of the same type. Note that 
this model still allows non-configurational changes where the number of type 
A interface sites and the number of type B interface sites may both be con- 
served after a growth event. Let fii(A) denote the number of type A interface 
sites adjoining site % and let fii(B) denote the number of type B interface sites 
adjoining site i. The transition probabilities for site i to change to an A site 
or a B site are 



Pi-, A' 



Pi 



kMa fii(A) 1 



kN a +Mb Ma 



n, 



n,(B) 1 



kMa+Mb M b n c 



(49) 
(50) 
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4-1 Mean Field Equilibrium 



Substituting (49,50) into the equilibrium condition (10) yields 

N N 

Kj2ni(A)6 i ,B = Y, fi i( B ) S i,A- (51) 



We can simplify this expression by noting that all type A sites adjoining a 
type B site are interface sites and vice-versa hence we have the identities 



N N 

Y / n i (A)5 i , B =J2 n i( A ) 5 i,B, (52) 

i i 
N N 

Y j n l {B)5 hA = Y j n t {B)5 lA . (53) 



Substituting the above identities together with the identity (17) we obtain the 
equilibrium condition 

k = 1. (54) 



For other values of k it is expected that the system will evolve until all sites 
are the same type; type A for k > 1 and type B for k < 1. 



The mean field equilibrium results have been found to provide a good de- 
scription of the long time populations in stochastic simulations of this model 
carried out on [0, 100] x [—100, 100] square lattice sites with the top half of 
the region initially comprised of abnormal cells and the bottom half initially 
comprised of normal cells. In Fig. 4 we show four different time snapshots for 
the case k — 2. 



As in the WB model the interface is highly irregular and the advantaged 
species completely dominates. However, the time taken for the extinction of 
the disadvantaged species is considerably less than for the WB model. The 
extinction times for the interface model over a range of k are shown in Fig. 
5. As in the original WB model the extinction time rapidly increases with 
decreasing \k — 1| and there is again a clear possibility for a steady state with 
both species at k — 1. Comparing Fig. 2 and Fig. 5 we see that for 10 4 cells 
the extinction time is an order of magnitude less in the interface model than 
in the WB model. 
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4-2 Regression 



In the case where initially we have a single type A site surrounded by type B 
sites the transition probabilities are 



PB - A = —A' (55) 

PA - B= ^h> (56) 

Pa^a = 0, (57) 

* M =db (58) 

The probability for regression without growth after m steps is 



^(m) = £(-— r Tmn ' ( 59 ) 



1 - 



J\k + 4/ (« + 4)' 

(60) 



3 

K+4 



K + 1 

The probability for regression without growth after infinitely many steps is 
again 

Pa^b(°o) = ——. (61) 

a T 1 

4-3 Limit k — > oo 



The k — > oo limit of this interface model is equivalent to the k — > oo limit of 
the WB model and is thus equivalent to the Eden B model. 



5 Selective Interface Model 



In the selective interface model an interface site is selected at random with a 
bias for selecting a type A site. Then one of the interface neighbouring sites 
of the selected site is chosen at random and converted to the same type as the 
initially selected site. The transition probabilities are 
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Mb 1 1^ 5i,, B 
P ^ B kMa + Mb k + 1 Mb V r* (A) ' [M) 

where the sum YTf denotes a sum over the neighbours of site %. 
5. 1 Mean Field Equilibrium 

The mean field equilibrium condition (10) for this model is 

a c c n c c 

r-» Oi',A r _V^V^ °V,B 



«EE^ = EE^. (64) 



This expression can be simplified. First we reverse the order of the sum over 
interface sites, an d the sum over neighbours, Y?n to obtain 

it ' nm ^(^) 



It is now a trivial matter to perform the sum over neighbours; 

§>, B = n*(£), (66) 

i' 

X;^,A = ni(A), (67) 

v 

leading to the equilibrium condition 

K M A = M B . (68) 



In general we might expect that the number of type A interface sites and 
the number of type B interface sites are approximately equal. The mean field 
steady state equilibrium that is consistent with this expectation is; k — 1 and 
Mb = Ma, or Mb = Ma = 0. The result Mb = Ma = is consistent with all 
type A sites for k > 1 and all type B sites for k < 1. 



5.2 Regression 



In the case where a single type A site is surrounded by n c type B sites the 
transition probabilities reduce to 
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(70) 



so that the probability for regression without growth in this model is 



4 

K+4' 



5.3 Limit k — > oo 



In the limit k — > oo the transition probability is 



1 ^ 6? A 



N A^-f n v {B) 



which is equivalent to the Eden C model [15]. 



6 Cross-Feeding Model 



Consider a colony of Ma type A interface cells and Af B type B interface cells. 
Select a cell at random with a bias for choosing type B. Convert the chosen 
cell to the opposite type. The transition probabilities are 



kNa 1 , f70 . 
P ^ A= .Af A+ Af B AT B 5 ^ (72) 

^ A) - B = KN A +N B W A kA - (73) 



Hence 



^ b = ^Ta^' (75) 

so that the probability for regression without growth is k jJ^1j^ b and the equi- 
librium condition is kMa = Nb- 

In the k — > oo limit the probability for a particular B cell to become an A cell 
is 1/Mb which is equivalent to the Eden A model [15]. 
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7 Cross-Feeding Model with Neighbour Bias 



Consider a colony of Ma type A interface cells and Mb type B interface cells. 
Select either the set of A cells or the set of B cells with a bias for choosing 
the set of B cells. Pick a cell at random from the selected set with a weighting 
according to the number of neighbouring cells of opposite type. Convert the 
chosen cell to the opposite type. The transition probabilities are 



kM a rii(A) 

^^^MaTm-b^ 6 ^ (76) 

P ^ b= kMa+M b M a &i * (77) 

where the B denotes an interface B site and the A denotes an interface A site. 
In this model the transition probabilities are not automatically normalized. 
However, it is a simple matter to calculate the normalization 



kMa Mb" N 



1 /kM a Mi 



+ ?]E».(%. (79) 



kM a +M b \M b M a 
where we have used the identity 



N 



N N 



J2n t (A)5 hB = J2n l (B)5 hA , (80) 



which follows immediately from the identity (17). 

Summing over all sites we obtain the mean field probability that a B cell 
becomes an A cell and vice-versa: 



K Af B + Ma 
A/s ^ Ma 

(83) 

The mean field equilibrium condition in this model is thus kMa 2 = Mb 2 - 
Again if we use the approximation Ma ~ Mb then the equilibrium condition 
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is satisfied by; either k — 1 and Ma ~ Ab, or A/a = A/b = 0. The latter case 
occurring when one or other species completely dominates. 

In this model if we start with a single type A cell surrounded by n c type B 
cells it is a simple matter to show that the probability for regression without 
growth is 

The k — > oo limit of this model reduces to the original Eden model [5,6], as 
in (48) above. 



8 Discussion 

In this paper we have introduced a class of stochastic growth models involving 
competition between two phases in which one of the phases has a competi- 
tive advantage k. This class of models includes the WB model for tumour 
growth. In the limit k — > oo where only one of the phases is unstable the 
class of models reduces to well known variants of the Eden model. We have 
derived mean field equilibrium results and regression probabilities (expressing 
the probability that the advantaged phase dies out) for the class of competitive 
growth models. The results are found to be in good agreement with stochastic 
simulations carried out on a square lattice with periodic boundaries. 

In the class of models studied here an equilibrium configuration with both 
types of cells can only occur when k — 1. An oft quoted result for the WB 
model is that the regression probability is ^ (or ^-)- Our simulations and 
analysis show that the result applies both to interface implementations of 
the WB model and to the full WB model. 

In subsequent work we plan to explore the important question of the rate of 
approach to the equilibrium state. 

This work has been supported by the Australian Research Council. 
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Fig. 1. Configurations of normal cells (white) and abnormal cells (black) in the WB 
model with carcinogenic advantage k = 2; (a) initially, (b) after 8 x 10 5 steps, (c) 
after 16 x 10 5 steps, (d) after 24 x 10 5 steps. 
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Fig. 2. Time taken for the extinction of normal cells (n > 1) and the extinction of 
abnormal cells (k < 1) as a function of logK in the WB model. 
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Fig. 3. The frequency of regression within (a) three (b) five and (c) ten steps versus 
the total number of runs used in the frequency estimate for the WB model. The 
lower horizontal line is the probability of regression after three steps and the two 
upper horizontal lines are the probabilities for regression without growth after (b) 
five and (c) ten steps. 
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Fig. 4. Configurations of normal cells (white) and abnormal cells (black) in the 
Interface Model model with carcinogenic advantage k = 2; (a) initially, (b) after 
25 x 10 3 steps, (c) after 50 x 10 3 steps, (d) after 75 x 10 3 steps. 
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Fig. 5. Time taken for the extinction of normal cells (n > 1) and the extinction of 
abnormal cells (k < 1) as a function of logK in the Interface Model. 
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